index_upload

1. Data preparation

library(dplyr)
library(traitstrap)
library(ggplot2)

morphometric_data <- utils::read.csv(here::here("data", "morphometric_data.csv"), sep = ";", header = T, dec = ".")

morpho_data <- morphometric_data %>%
  select(-c(variable_abbreviation, variable_unit)) %>%
  t() %>%
  as.data.frame() %>%
  janitor::row_to_names(row_number = 1) %>%
  `rownames<-`(NULL)%>%
  # delete for now (n=1) and Eurypharynx pelecanoides present 7 missing traits 
  filter(!species%in% c("Diaphus_sp","Eurypharynx_pelecanoides"))

# replace empty value by NA 
morpho_data[morpho_data ==""] <- NA

# Numeric variables
morpho_data[, 4:23] <- sapply(morpho_data[, 4:23], as.numeric)

1.1.Data imputation

• mice algorithm: n imputation = 5, n iterations = 50

#select numeric variables for imputation 
original_data <- morpho_data %>%
   filter(species!="Eurypharynx_pelecanoides") %>% 
  select(1:23)

imputation <-
  mice::mice(
    original_data,
    m = 5,
    maxit = 50,
    printFlag = F
  )

imputed_data <- mice::complete(imputation)

1.2. Species * traits

• calculate functional traits

# calculate functional numeric traits
numeric_traits <- imputed_data %>%
  na.omit() %>%
  select(-individual_code) %>%
  mutate(
    eye_size = eye_diameter / head_depth,
    orbital_length = eye_diameter / standard_length,
    oral_gape_surface = mouth_width * mouth_depth / body_width * body_depth,
    oral_gape_shape = mouth_depth / mouth_width,
    oral_gape_position = distance_upper_jaws_bottom_head / head_depth,
    lower_jaw_length = lower_jaw_length / standard_length,
    head_length = head_length / standard_length,
    body_depth = body_depth / standard_length,
    pectoral_fin_position = distance_pectoral_bottom_body / body_depth_pectoral_insertion,
    pectoral_fin_insertion = prepectoral_length / standard_length,
    transversal_shape = body_depth / body_width,
    dorsal_fin_insertion = predorsal_length / standard_length,
    eye_position = eye_height / head_depth,
    operculum_volume = operculum_depth / operculum_width,
    gill_outflow = operculum_width,
    caudal_throttle_width = caudal_peduncle_min_depth
  ) %>%
  select(
    species,
    eye_size,
    orbital_length,
    gill_outflow,
    oral_gape_surface,
    oral_gape_shape,
    oral_gape_position,
    lower_jaw_length,
    head_length,
    body_depth,
    pectoral_fin_position,
    pectoral_fin_insertion,
    transversal_shape,
    caudal_throttle_width,
    dorsal_fin_insertion,
    eye_position,
    operculum_volume
  ) %>%
  group_by(species) %>%
  summarise(across(everything(), mean, na.rm = TRUE)) %>%
  arrange(species)

# categorical traits for species without NA
cat_morpho <- morpho_data %>%
  select(
    species,
    ventral_photophores,
    gland_head,
    chin_barbel,
    small_teeth,
    large_teeth,
    fang_teeth,
    retractable_teeth,
    internal_teeth,
    gill_raker_types,
    oral_gape_axis
  ) %>%
    na.omit() %>%
  distinct() %>%
  arrange(species)

# combined the two data frames
fish_traits <- numeric_traits %>%
  inner_join(cat_morpho, by = "species") %>%
  arrange(species) %>% 
  mutate(species = stringr::str_replace(species,  "Cyclothone_sp","Cyclothone")) %>% 
  filter(species!="Eurypharynx_pelecanoides") %>% 
  tibble::column_to_rownames("species")%>%
  # assign trait type 
  # as.factor for qualitative traits
  mutate_if(is.character, as.factor)%>%
  # as.ordered for ordered variables
  mutate_at(c("gill_raker_types", "oral_gape_axis"), as.ordered)

Data summary

fish_traits_sum <- fish_traits %>% 
  mutate(across(where(is.numeric), round, 2))

## Display the table ----
htmltools::tagList(DT::datatable(fish_traits_sum))

traits correlation

M <-cor(numeric_traits[, c(-1)])

ggcorrplot::ggcorrplot(M, hc.order = TRUE, type = "lower",
                       lab = TRUE, tl.cex = 9, lab_size = 3)

# list of species 
sp_names <- c(rownames(fish_traits))

# taxonomic_families
taxonomic_families <- sp_names %>%
  as.data.frame() %>%
  `colnames<-`("species") %>% 
  mutate(
    family = case_when(
      species %in%
        c(
          "Benthosema_glaciale",
          "Ceratoscopelus_maderensis",
          "Diaphus_metopoclampus",
          "Lampanyctus_ater",
          "Lampanyctus_crocodilus",
          "Lampanyctus_macdonaldi",
          "Lobianchia_gemellarii",
          "Myctophum_punctatum",
          "Notoscopelus_bolini",
          "Notoscopelus_kroyeri",
          "Bolinichthys_supralateralis"
        ) ~ "Myctophidae",
      species %in% c(
        "Borostomias_antarcticus",
        "Chauliodus_sloani",
        "Malacosteus_niger",
        "Melanostomias_bartonbeani",
        "Stomias_boa"
      ) ~ "Stomiidae",
      species %in% c(
        "Holtbyrnia_anomala",
        "Holtbyrnia_macrops",
        "Maulisia_argipalla",
        "Maulisia_mauli",
        "Maulisia_microlepis",
        "Normichthys_operosus",
        "Searsia_koefoedi",
        "Sagamichthys_schnakenbecki"
      ) ~ "Platytroctidae",
      species %in% c("Sigmops_bathyphilus",
                     "Gonostoma_elongatum") ~ "Gonostomatidae",
      species %in% c(
        "Argyropelecus_hemigymnus",
        "Maurolicus_muelleri",
        "Argyropelecus_olfersii"
      ) ~ "Sternoptychidae",
      species == "Anoplogaster_cornuta" ~ "Anoplogastridae",
      species %in% c("Arctozenus_risso", "Paralepis_coregonoides") ~ "Paralepididae",
      species == "Bathylagus_euryops" ~ "Bathylagidae",
      species == "Cyclothone" ~ "Gonostomatidae",
      species == "Derichthys_serpentinus" ~ "Derichthyidae",
      species == "Eurypharynx_pelecanoides" ~ "Eurypharyngidae",
      species == "Evermannella_balbo" ~ "Evermannellidae",
      species == "Lestidiops_sphyrenoides" ~ "Lestidiidae",
      species == "Melanostigma_atlanticum" ~ "Zoarcidae",
      species %in% c("Photostylus_pycnopterus",
                     "Xenodermichthys_copei") ~ "Alepocephalidae",
      species == "Serrivomer_beanii" ~ "Serrivomeridae"
    )
  )%>% 
  mutate(
    order = case_when(
      family =="Myctophidae" ~ "Myctophiformes",
      family %in% c("Stomiidae","Gonostomatidae", "Sternoptychidae") ~  "Stomiiformes",
      family %in% c("Platytroctidae","Alepocephalidae") ~ "Alepocephaliformes",
      family == "Anoplogastridae" ~ "Trachichthyiformes",
      family %in% c("Paralepididae","Evermannellidae","Lestidiidae") ~ "Aulopiformes",
      family ==  "Bathylagidae" ~ "Argentiniformes",
      family %in% c("Derichthyidae","Serrivomeridae") ~ "Anguilliformes",
      family ==  "Eurypharyngidae" ~"Saccopharyngiformes",
      family == "Zoarcidae" ~ "Perciformes",
    )
  )

1.3. Species * assemblages matrix

Number of trawl hauls per depth

• Epipelagic = 8
• Upper mesopelagic = 26
• Lower mesopelagic = 16
• Bathypelagic = 16

# Metadata
metadata <-  utils::read.csv(here::here("data", "metadata.csv"), sep = ";", header = T, dec = ".")%>%
  # calculation of standardized biomass values (vertical  trawl opening * horizontal trawl opening * distance traveled)  
  mutate(volume_filtered = 24*58*distance_trawled_m)

ggplot(metadata, aes(x=depth))  +
  ylab ("Number of trawls")+
  xlab("Immersion depth (m)")+
  geom_histogram(binwidth=100, col="white", fill=alpha("black",0.55))+
  theme_light()+
  coord_flip()+ 
  scale_x_reverse()+
  labs(fill= "")+
  guides(fill="none")+
  scale_y_continuous(breaks = c(2,4,6,8,10,12))+
  theme(axis.text.x= element_text(size=12),
        axis.text.y= element_text(size=12),
        axis.title.y = element_text( size=12),
        axis.ticks = element_blank())

Biomass data

biomass_sum <- depth_fish_biomass %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "depth_layer") %>%
  tidyr::pivot_longer(!depth_layer, names_to = "species", values_to = "biomass") %>%
  group_by(depth_layer) %>%
  filter(biomass > 0) %>%
  summarise(biomass_depth = round(sum(biomass), 3),
            n = n())

htmltools::tagList(DT::datatable(biomass_sum))

1.4. Traits types

The first column contains traits name. The second column contains traits type following this code:

• N: nominal trait (factor variable)
• O: ordinal traits (ordered variable)
• Q: quantitative traits (numeric values)

fish_traits_cat <- utils::read.csv(here::here("data", "fish_traits_cat.csv"), sep = ";", header = T, dec = ".") 
htmltools::tagList(DT::datatable(fish_traits_cat))

2.Functional spaces

2.1 Compute data summaries

## Summary of the assemblages * species data.frame ----
asb_sp_fish_summ <- mFD::asb.sp.summary(asb_sp_w = depth_fish_biomass)
asb_sp_fish_occ  <- asb_sp_fish_summ$"asb_sp_occ"

htmltools::tagList(DT::datatable(asb_sp_fish_occ))

2.2 Computing distances between species based on functional traits

• We have non-continuous traits so we use the Gower distance (metric = “gower”) as this method allows traits weighting.
• scale_euclid = TRUE

sp_dist_fish <- mFD::funct.dist(
  sp_tr         = fish_traits,
  tr_cat        = fish_traits_cat,
  metric        = "gower",
  scale_euclid  = "scale_center",
  ordinal_var   = "classic",
  weight_type   = "equal",
  stop_if_NA    = TRUE)

2.3 Building functional spaces and chosing the best one

2.3.1 Computing several multimensional functional spaces and assessing their quality

• mFD evaluates the quality of PCoA-based multidimensional spaces according to the deviation between trait-based distances and distances in the functional space (extension of Maire et al. (2015) framework).

fspaces_quality_fish <- mFD::quality.fspaces(
  sp_dist             = sp_dist_fish,
  maxdim_pcoa         = 10,
  deviation_weighting = "absolute",
  fdist_scaling       = FALSE,
  fdendro             = "average")

## Quality metrics of functional spaces ----
round(fspaces_quality_fish$"quality_fspaces", 3)

               mad
pcoa_1d      0.140
pcoa_2d      0.077
pcoa_3d      0.048
pcoa_4d      0.029
pcoa_5d      0.022
pcoa_6d      0.016
pcoa_7d      0.015
pcoa_8d      0.015
pcoa_9d      0.017
pcoa_10d     0.020
tree_average 0.040

Variance explained by each axis

# Extract eigenvalues information
eigenvalues_info <- fspaces_quality_fish$"details_fspaces"$"pc_eigenvalues"

# Create a dataframe to store the results
variance_df <- data.frame(
  PC = c("PC1", "PC2", "PC3", "PC4"),
  VarianceExplained = c(
    eigenvalues_info[1, "Cum_corr_eig"] * 100,
    (eigenvalues_info[2, "Cum_corr_eig"] - eigenvalues_info[1, "Cum_corr_eig"]) * 100,
    (eigenvalues_info[3, "Cum_corr_eig"] - eigenvalues_info[2, "Cum_corr_eig"]) * 100,
    (eigenvalues_info[4, "Cum_corr_eig"] - eigenvalues_info[3, "Cum_corr_eig"]) * 100
  )
) %>% 
  mutate(across(where(is.numeric), round, 2))

htmltools::tagList(DT::datatable(variance_df))

2.3.2 Illustrating the quality of the functional spaces

mFD::quality.fspaces.plot(
  fspaces_quality            = fspaces_quality_fish,
  quality_metric             = "mad",
  fspaces_plot               = c("tree_average", "pcoa_2d", "pcoa_3d", 
                                 "pcoa_4d", "pcoa_5d", "pcoa_6d"),
  name_file                  = NULL,
  range_dist                 = NULL,
  range_dev                  = NULL,
  range_qdev                 = NULL,
  gradient_deviation         = c(neg = "darkblue", nul = "grey80", pos = "darkred"),
  gradient_deviation_quality = c(low = "yellow", high = "red"),
  x_lab                      = "Trait-based distance")

2.3.3 Testing the correlation between functional axes and traits

sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"

# As we have 26 traits we have to split the df to see correlation between functional axes and traits 
# first set ----
fish_traits_1 <- fish_traits%>%
  select(1:9)

fish_tr_faxes <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits_1, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = T)

## Print traits with significant effect ----
fish_tr_faxes$"tr_faxes_stat"[which(fish_tr_faxes$"tr_faxes_stat"$"p.value" < 0.05), ]

                trait axis         test stat value p.value
1            eye_size  PC1 Linear Model   r2 0.161  0.0093
2            eye_size  PC2 Linear Model   r2 0.183  0.0052
3            eye_size  PC3 Linear Model   r2 0.133  0.0192
4            eye_size  PC4 Linear Model   r2 0.139  0.0163
5      orbital_length  PC1 Linear Model   r2 0.410  0.0000
10       gill_outflow  PC2 Linear Model   r2 0.511  0.0000
13  oral_gape_surface  PC1 Linear Model   r2 0.139  0.0165
14  oral_gape_surface  PC2 Linear Model   r2 0.513  0.0000
18    oral_gape_shape  PC2 Linear Model   r2 0.166  0.0083
24 oral_gape_position  PC4 Linear Model   r2 0.221  0.0019
26   lower_jaw_length  PC2 Linear Model   r2 0.699  0.0000
29        head_length  PC1 Linear Model   r2 0.325  0.0001
30        head_length  PC2 Linear Model   r2 0.376  0.0000
34         body_depth  PC2 Linear Model   r2 0.352  0.0000
35         body_depth  PC3 Linear Model   r2 0.203  0.0031

## Plot ----
fish_tr_faxes$"tr_faxes_plot"

# second set ----
fish_traits_2 <- fish_traits%>%
  select(10:18)

fish_tr_faxes_2 <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits_2, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = T)

## Print traits with significant effect ----
fish_tr_faxes_2$"tr_faxes_stat"[which(fish_tr_faxes_2$"tr_faxes_stat"$"p.value" < 0.05), ]

                    trait axis           test stat value p.value
2   pectoral_fin_position  PC2   Linear Model   r2 0.175  0.0065
5  pectoral_fin_insertion  PC1   Linear Model   r2 0.279  0.0004
6  pectoral_fin_insertion  PC2   Linear Model   r2 0.388  0.0000
9       transversal_shape  PC1   Linear Model   r2 0.123  0.0246
11      transversal_shape  PC3   Linear Model   r2 0.220  0.0020
14  caudal_throttle_width  PC2   Linear Model   r2 0.410  0.0000
19   dorsal_fin_insertion  PC3   Linear Model   r2 0.184  0.0051
23           eye_position  PC3   Linear Model   r2 0.204  0.0030
32    ventral_photophores  PC4 Kruskal-Wallis eta2 0.590  0.0000
33             gland_head  PC1 Kruskal-Wallis eta2 0.245  0.0011

## Plot ----
fish_tr_faxes_2$"tr_faxes_plot"

# third set ----
fish_traits_3 <- fish_traits%>%
  select(19:26)

fish_tr_faxes_3 <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits_3, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = T)

## Print traits with significant effect ----
fish_tr_faxes_3$"tr_faxes_stat"[which(fish_tr_faxes_3$"tr_faxes_stat"$"p.value" < 0.05), ]

              trait axis           test stat value p.value
1       chin_barbel  PC1 Kruskal-Wallis eta2 0.142  0.0107
4       chin_barbel  PC4 Kruskal-Wallis eta2 0.183  0.0043
5       small_teeth  PC1 Kruskal-Wallis eta2 0.317  0.0003
9       large_teeth  PC1 Kruskal-Wallis eta2 0.461  0.0000
10      large_teeth  PC2 Kruskal-Wallis eta2 0.205  0.0027
13       fang_teeth  PC1 Kruskal-Wallis eta2 0.410  0.0000
21   internal_teeth  PC1 Kruskal-Wallis eta2 0.082  0.0408
23   internal_teeth  PC3 Kruskal-Wallis eta2 0.619  0.0000
25 gill_raker_types  PC1 Kruskal-Wallis eta2 0.329  0.0007
26 gill_raker_types  PC2 Kruskal-Wallis eta2 0.381  0.0003
30   oral_gape_axis  PC2 Kruskal-Wallis eta2 0.360  0.0004
31   oral_gape_axis  PC3 Kruskal-Wallis eta2 0.240  0.0038

## Plot ----
fish_tr_faxes_3$"tr_faxes_plot"

Summary of traits with a significant effect

sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"

fish_tr_faxes <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = F)

## Print traits with significant effect ----
traits_effect <- fish_tr_faxes[which(fish_tr_faxes$p.value< 0.05),] %>% 
  as.data.frame() %>% 
  arrange(axis, desc(value))

htmltools::tagList(DT::datatable(traits_effect))

2.4 Plotting the selected functional space and position of species

# Convert sp_faxes_coord_fish to a data frame
sp_coord_community <- as.data.frame(sp_faxes_coord_fish[, c("PC1", "PC2")]) %>%
  tibble::rownames_to_column(var = "species")

# Depth_layer ----
# Transform biomass data
sp_all_layers <- depth_fish_biomass %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "depth_layer") %>%
  tidyr::pivot_longer(!depth_layer, values_to = "biomass", names_to = "species") %>%
  filter(biomass > 0)

# Create a list to store data frames for each depth layer
all_data_layers <- lapply(unique(sp_all_layers$depth_layer), function(layer) {
  sp_layer <- sp_all_layers %>%
    filter(depth_layer == layer)
  
  # Mark presence or absence
  sp_coord_layer <- sp_coord_community %>%
    mutate(layer_presence = ifelse(species %in% sp_layer$species, "yes", "no")) %>%
    mutate(depth_layer = layer)
  
  return(sp_coord_layer)
})

# Combine all layers into one data frame
sp_coord_community_all_layers <- bind_rows(all_data_layers)

# Function to calculate the convex hull for a given data frame
calculate_hull <- function(data) {
  if (nrow(data) < 3) return(data)  # Handle cases with fewer than 3 points
  data %>%
    slice(chull(PC1, PC2))
}

# Initialize lists to store hulls
hull_all_combined <- NULL
hull_layer_combined <- NULL

for(layer in unique(sp_coord_community_all_layers$depth_layer)) {
  sp_layer <- sp_coord_community_all_layers %>%
    filter(depth_layer == layer)
  
  hull_all <- calculate_hull(sp_layer) %>%
    mutate(depth_layer = layer, type = "all")
  
  hull_layer <- calculate_hull(sp_layer %>% filter(layer_presence == "yes")) %>%
    mutate(depth_layer = layer, type = "present")
  
  hull_all_combined <- bind_rows(hull_all_combined, hull_all)
  hull_layer_combined <- bind_rows(hull_layer_combined, hull_layer)
}

# Define depth layer levels in desired order
depth_levels <- c("Epipelagic", "Upper mesopelagic", "Lower mesopelagic", "Bathypelagic")

# Ensure depth_layer is a factor with specified levels
sp_coord_community_all_layers$depth_layer <- factor(sp_coord_community_all_layers$depth_layer, levels = depth_levels)
hull_all_combined$depth_layer <- factor(hull_all_combined$depth_layer, levels = depth_levels)
hull_layer_combined$depth_layer <- factor(hull_layer_combined$depth_layer, levels = depth_levels)

# Define colors for each depth layer
depth_colors <- c("Epipelagic" = "#FEA520", 
                  "Upper mesopelagic" = "#D62246", 
                  "Lower mesopelagic" = "#6255B4", 
                  "Bathypelagic" = "#3C685A")

# Plotting
plot_depth <- ggplot() +
  # Plot present species
  geom_point(data = sp_coord_community_all_layers %>% filter(layer_presence == "yes"), 
             aes(PC1, PC2, col = depth_layer), size = 2, shape=19) +
  # Plot absent species
  geom_point(data = sp_coord_community_all_layers %>% filter(layer_presence == "no"), 
             aes(PC1, PC2), col = "grey", shape = 8, size = 2) +
  # Plot hulls
  geom_polygon(data = hull_all_combined, aes(x = PC1, y = PC2), 
               fill =  "gray70", alpha = 0.2) +
  geom_polygon(data = hull_layer_combined, aes(x = PC1, y = PC2, group = depth_layer, fill = depth_layer), 
               alpha = 0.4) +
  theme_light() +
  scale_color_manual(values = depth_colors) +
  scale_fill_manual(values = depth_colors) +
  guides(col = "none", shape = "none", fill = "none") +
  facet_wrap(~ depth_layer, nrow=1 ) +  
  theme(
    strip.text.x = element_text(size = 14, color = "gray20"),
    strip.background = element_rect(fill = "white"),
    aspect.ratio = 1,
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank()
  )


beta_fd_indices_fish <- mFD::beta.fd.multidim(
  sp_faxes_coord   = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")],
  asb_sp_occ       = asb_sp_fish_occ,
  check_input      = TRUE,
  beta_family      = c("Jaccard"),
  details_returned = T)

Serial computing of convex hulls shaping assemblages with conv1

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
Serial computing of intersections between pairs of assemblages with inter_geom_coord

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |======================================================================| 100%

vertices_community <-beta_fd_indices_fish$"details"$"asb_vertices"$Bathypelagic 

sp_coord_community <- as.data.frame(sp_faxes_coord_fish[, 1:2]) %>% 
  tibble::rownames_to_column(var = "species") %>% 
  mutate(vertices=case_when(species%in%vertices_community~"vertice",
                            !species%in%vertices_community~"not_vertice")) 

hull <- sp_coord_community%>% 
  slice(chull(PC1, PC2))

plot_com <- ggplot(data = sp_coord_community,aes(PC1, PC2)) +
  scale_color_manual(values = c("grey", "black"))+
  geom_point(size = 2.5, alpha=0.8, aes( shape=vertices, col= vertices)) +
  geom_polygon(data=hull, aes(x = PC1, y = PC2),
               fill = "gray70", alpha = 0.2) +
  theme_light() +
  geom_hline(yintercept = 0, linetype = 2, col = "gray45") +
  geom_vline(xintercept = 0, linetype = 2, col = "gray45") +
  guides(col = "none", shape = "none", fill = "none") +
  theme(
    strip.text.x = element_text(size = 14, color = "gray20"),
    strip.background = element_rect(fill = "white"),
    aspect.ratio = 1,
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank()
  )

ggpubr::ggarrange(
  labels = c("A","B"),
  plot_com,   
  plot_depth, 
  nrow = 2
) 

ggsave("functional_space.png", path = "figures", dpi = 800, height = 8, width = 10)

3. SES Functional diversity indices

3.1.calculate and plot SES

• Standard Effect Size (SES): to eliminate the influence of species richness on the functional diversity indices (Mouchet et al., 2010). Measures the deviation from the random expectation in standard deviation units

• null model frequency: Randomize community data matrix abundances (here biomasss) within species (maintains species occurrence frequency)

• FRic Functional Richness: the proportion of functional space filled by species of the studied assemblage, i.e. the volume inside the convex-hull shaping species. To compute FRic the number of species must be at least higher than the number of functional axis + 1.

• FDis Functional Dispersion: the biomass weighted deviation of species traits values from the center of the functional space filled by the assemblage i.e. the biomass-weighted mean distance to the biomass-weighted mean trait values of the assemblage.

• FDiv Functional Divergence: the proportion of the biomass supported by the species with the most extreme functional traits i.e. the ones located close to the edge of the convex-hull filled by the assemblage.

• FEve Functional Evenness: the regularity of biomass distribution in the functional space using the Minimum Spanning Tree linking all species present in the assemblage.

• calcul des indices avec package FD

• correction distance matrix (obtenue avec Gower) = “none” or “lingoes” transformation avec “sqrt” ne permet pas de transformer la matrice en “euclidean”

• randomisation des traits des espèces par couche de profondeur

# # Model parameters ----
# # Number of simulations
# n_simulations <- 9
# 
# # correction method to use when the distance matrix cannot be represented in a Euclidean space
# corr_method <- "lingoes"# "none"
# 
# # Depth layers
# depth_layers <- rownames(depth_fish_biomass)
# 
# # List of indices to be calculated
# indices <- c("FRic", "FDis", "FDiv", "FEve")
# 
# # Matrices for storing observed and simulated results
# dbFD_result_obs <- matrix(NA, nrow = length(depth_layers), ncol = length(indices),
#                           dimnames = list(depth_layers, indices))
# 
# dbFD_result_sim <- array(NA, dim = c(length(depth_layers), n_simulations, length(indices)),
#                          dimnames = list(depth_layers, paste0("Sim.", 1:n_simulations), indices))
# 
# # Function to randomize traits while preserving factor levels----
# randomize_traits <- function(traits_mat) {
#   randomized_mat <- traits_mat
#   # Randomisation of columns
#   for (trait in colnames(traits_mat)) {
#     randomized_mat[, trait] <- sample(traits_mat[, trait], replace = FALSE)
#   }
#   return(randomized_mat)
# }
# 
# # Loop on each depth layer ----
# for (layer in depth_layers) {
# 
#   # Select the species present in the selected depth layer
#   species_in_layer <- colnames(depth_fish_biomass)[depth_fish_biomass[layer, ] > 0]
# 
#   # Filter lines (traits) and remove constant columns
#   traits_layer <- fish_traits[species_in_layer, , drop = FALSE] %>%
#     select(where(~ length(unique(.)) > 1))
# 
#   # If no variable lines remain, we move on to the next layer
#   if (ncol(traits_layer) == 0) next
# 
#   # Correct formatting of biomass_layer with depth layer in rownames
#   # Converts biomass data to matrix, which is required by FD::dbFD()
#   biomass_layer <- matrix(as.numeric(depth_fish_biomass[layer, species_in_layer, drop = FALSE]),
#                           nrow = 1, dimnames = list(layer, species_in_layer))
# 
#   # Calculation of observed indices with FD::dbFD
#   dbFD_result <- FD::dbFD(x = traits_layer, w.abun = TRUE, m = 4, a = biomass_layer, messages=F,
#                           corr = corr_method)
# 
#   # Storage of observed values
#   dbFD_result_obs[layer, "FRic"] <- dbFD_result$FRic
#   dbFD_result_obs[layer, "FDis"] <- dbFD_result$FDis
#   dbFD_result_obs[layer, "FEve"] <- dbFD_result$FEve
#   dbFD_result_obs[layer, "FDiv"] <- dbFD_result$FDiv
# 
#   # Simulations for each layer
#   for (sim in 1:n_simulations) {
#     randomized_traits <- randomize_traits(traits_layer)
# 
#     dbFD_sim_result <- FD::dbFD(x = randomized_traits, w.abun = TRUE, m = 4, a = biomass_layer, messages=F,
#                                 corr = corr_method)
# 
#     # Storage of simulated values
#     dbFD_result_sim[layer, sim, "FRic"] <- dbFD_sim_result$FRic
#     dbFD_result_sim[layer, sim, "FDis"] <- dbFD_sim_result$FDis
#     dbFD_result_sim[layer, sim, "FEve"] <- dbFD_sim_result$FEve
#     dbFD_result_sim[layer, sim, "FDiv"] <- dbFD_sim_result$FDiv
#   }
# }
# 
# # Calculation of the SES (Standardized Effect Size) ----
# SES_results <- data.frame(depth_layer = rownames(dbFD_result_obs))
# 
# for (index in indices) {
#   meanNullFD <- rowMeans(dbFD_result_sim[, , index], na.rm = TRUE)
#   sdNullFD <- apply(dbFD_result_sim[, , index], 1, sd, na.rm = TRUE)
#   SES_results[[index]] <- (dbFD_result_obs[, index] - meanNullFD) / sdNullFD
# }
# 
# # Plot ----
# results_df <- tidyr::pivot_longer(
#   SES_results,
#   cols = -depth_layer,
#   names_to = "index",
#   values_to = "SES_FD"
# )
# 
# results_df$depth_layer <- factor(
#   results_df$depth_layer,
#   levels = c(
#     "Epipelagic",
#     "Upper mesopelagic",
#     "Lower mesopelagic",
#     "Bathypelagic"
#   )
# )
# 
# results_df$index <- factor(
#   results_df$index,
#   levels = c("FRic",
#              "FDis",
#              "FDiv",
#              "FEve"),
#   labels = c(
#     "Functional richness",
#     "Functional dispersion",
#     "Functional divergence",
#     "Functional evenness"
#   )
# )
# 
# ggplot(results_df, aes(x = depth_layer, y = SES_FD, color = depth_layer)) +
#   geom_point(size = 3) +
#   facet_wrap(~index) +
#   geom_hline(yintercept = c(1.96, -1.96), linetype = "dashed", color = "gray40", size=0.8) +
#   scale_color_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
#   labs(x = "", y = "Standard Effect Size (SES)") +
#   theme_light() +
#   theme(axis.text.x = element_blank(),
#         strip.text.x = element_text(size = 14, color = "black"),
#         strip.background = element_rect(fill = "white"),
#         axis.title = element_text(size = 13),
#         axis.text = element_text(size = 13)) +
#   guides(col="none", fill="none")
# 
# ggsave("SES_dbFD.png", path = "figures", dpi = 700, height = 5, width = 7)

3.2. Two levels (stations and depth layer)

station_sp <- rbind(data_biomass_2002_2019, data_biomass_2021_2022) %>%
  as.data.frame() %>%
  left_join(metadata) %>%
  select(species, biomass_sp, volume_filtered, station) %>%
  # Divide biomass by the volume filtered at each trawl (g.m3)
  mutate(biomass_cpu = (biomass_sp / volume_filtered) * 1000) %>%
  select(species, biomass_cpu, station) %>%
  replace(is.na(.), 0) %>%
  group_by(species, station) %>%
  mutate(biomass = sum(biomass_cpu)) %>%
  select(-biomass_cpu) %>%
  distinct() %>%
  tidyr::pivot_wider(names_from = species, values_from = biomass) %>%
  replace(is.na(.), 0) %>%
  arrange(station) %>%
  filter(!station %in% c("H0411", "L0731", "L0736")) %>%
  tibble::column_to_rownames(var = "station") %>%
  select(order(colnames(.))) %>%
  as.matrix()


# Model parameters ----
n_simulations <- 9

# Correction method for non-Euclidean distances
corr_method <- "lingoes"
depth_layers <- rownames(station_sp)  

# Indices to calculate
indices <- c("FRic", "FDis", "FDiv", "FEve")

# Function to randomize traits while preserving factor levels----
randomize_traits <- function(traits_mat) {
  randomized_mat <- traits_mat
  # Randomisation of columns
  for (trait in colnames(traits_mat)) {
    randomized_mat[, trait] <- sample(traits_mat[, trait], replace = FALSE)
  }
  return(randomized_mat)
}

# Data definition: list with depth and station matrices ----
analysis_levels <- list(
  depth_layer = depth_fish_biomass,
  station = station_sp
)

# Store results ----
all_results <- list()
dbFD_result_sim_list <- list()

# Function to calculate FD indices and SES ----
calculate_FD_and_SES <- function(data_matrix, traits_data, n_simulations, indices, corr_method) {
  layers <- rownames(data_matrix)
  dbFD_result_obs <- matrix(NA, nrow = length(layers), ncol = length(indices), 
                            dimnames = list(layers, indices))
  dbFD_result_sim <- array(NA, dim = c(length(layers), n_simulations, length(indices)),
                           dimnames = list(layers, paste0("Sim.", 1:n_simulations), indices))
  
  # Loop through each layer
  for (layer in layers) {
    species_in_layer <- colnames(data_matrix)[data_matrix[layer, ] > 0]
    traits_layer <- traits_data[species_in_layer, , drop = FALSE] %>%
      select(where(~ length(unique(.)) > 1)) %>%
      droplevels()
    
    if (ncol(traits_layer) == 0) next
    
    biomass_layer <- matrix(as.numeric(data_matrix[layer, species_in_layer, drop = FALSE]),
                            nrow = 1, dimnames = list(layer, species_in_layer))
    
    # Calculate FD indices for observed data
    dbFD_result <- FD::dbFD(x = traits_layer, w.abun = TRUE, m = 4, a = biomass_layer,
                            messages = FALSE, corr = corr_method)
    
    # Verify structure of dbFD_result
    indices_names <- names(dbFD_result)  
    
    # Extract valid indices
    valid_indices <- intersect(indices, indices_names)
    dbFD_result_obs[layer, valid_indices] <- dbFD_result[valid_indices]
    
    # Random simulations
    for (sim in 1:n_simulations) {
      randomized_traits <- randomize_traits(traits_layer)
      dbFD_sim_result <- FD::dbFD(x = randomized_traits, w.abun = TRUE, m = 4, a = biomass_layer,
                                  messages = FALSE, corr = corr_method)
      
      # Verify structure of dbFD_sim_result
      indices_names_sim <- names(dbFD_sim_result)
      valid_indices_sim <- intersect(indices, indices_names_sim)
      
      # Store simulation results
      dbFD_result_sim[layer, sim, valid_indices_sim] <- dbFD_sim_result[valid_indices_sim]
    }
  }
  
  # Calculate SES (Standardized Effect Size) ----
  SES_results <- data.frame(name = rownames(dbFD_result_obs))
  for (index in indices) {
    meanNullFD <- rowMeans(dbFD_result_sim[, , index], na.rm = TRUE)
    sdNullFD <- apply(dbFD_result_sim[, , index], 1, sd, na.rm = TRUE)
    SES_results[[index]] <- (dbFD_result_obs[, index] - meanNullFD) / sdNullFD
  }
  
  return(list(dbFD_result_obs = dbFD_result_obs, dbFD_result_sim = dbFD_result_sim, SES_results = SES_results))
}

# Loop through analysis levels
for (level in names(analysis_levels)) {
  data_matrix <- analysis_levels[[level]]
  layers <- rownames(data_matrix)
  
  # Matrices for observed and simulated results
  dbFD_result_obs <- matrix(NA, nrow = length(layers), ncol = length(indices),
                            dimnames = list(layers, indices))
  
  dbFD_result_sim <- array(NA, dim = c(length(layers), n_simulations, length(indices)),
                           dimnames = list(layers, paste0("Sim.", 1:n_simulations), indices))
  
  # Loop through each layer
  for (layer in layers) {
    species_in_layer <- colnames(data_matrix)[data_matrix[layer, ] > 0]
    traits_layer <- fish_traits[species_in_layer, , drop = FALSE] %>%
      select(where(~ length(unique(.)) > 1)) %>%
      droplevels()
    
    if (ncol(traits_layer) == 0) next
    
    biomass_layer <- matrix(as.numeric(data_matrix[layer, species_in_layer, drop = FALSE]),
                            nrow = 1, dimnames = list(layer, species_in_layer))
    
    # Calculate FD indices for the observed data
    dbFD_result <- FD::dbFD(x = traits_layer, w.abun = TRUE, m = 4, a = biomass_layer,
                            messages = FALSE, corr = corr_method)
    
    # Store observed results
    dbFD_result_obs[layer, "FRic"] <- dbFD_result$FRic
    dbFD_result_obs[layer, "FDis"] <- dbFD_result$FDis
    dbFD_result_obs[layer, "FEve"] <- dbFD_result$FEve
    dbFD_result_obs[layer, "FDiv"] <- dbFD_result$FDiv
    
    # Perform random simulations
    for (sim in 1:n_simulations) {
      randomized_traits <- randomize_traits(traits_layer)
      
      dbFD_sim_result <- FD::dbFD(x = randomized_traits, w.abun = TRUE, m = 4, a = biomass_layer,
                                  messages = FALSE, corr = corr_method)
      
      # Store simulated results
      dbFD_result_sim[layer, sim, "FRic"] <- dbFD_sim_result$FRic
      dbFD_result_sim[layer, sim, "FDis"] <- dbFD_sim_result$FDis
      dbFD_result_sim[layer, sim, "FEve"] <- dbFD_sim_result$FEve
      dbFD_result_sim[layer, sim, "FDiv"] <- dbFD_sim_result$FDiv
    }
  }
  
  # Save simulated results for normality tests
  dbFD_result_sim_list[[level]] <- dbFD_result_sim
  
  # Calculate SES (Standardized Effect Size)
  SES_results <- data.frame(level = level, name = rownames(dbFD_result_obs))
  for (index in indices) {
    meanNullFD <- rowMeans(dbFD_result_sim[, , index], na.rm = TRUE)
    sdNullFD <- apply(dbFD_result_sim[, , index], 1, sd, na.rm = TRUE)
    SES_results[[index]] <- (dbFD_result_obs[, index] - meanNullFD) / sdNullFD
  }
  
  all_results[[level]] <- SES_results
}

# Normality and symmetry tests ----
library(nortest)
library(moments)

# Normality and symmetry results
normality_results <- list()

for (level in names(all_results)) {
  SES_data <- all_results[[level]]
  simulated_values <- dbFD_result_sim_list[[level]]
  
  normality_test <- data.frame(level = level, name = SES_data$name)
  
  for (index in indices) {
    # Normality test (Shapiro-Wilk)
    normality_test[[paste0(index, "_shapiro_p")]] <- apply(simulated_values[, , index], 1, function(x) shapiro.test(x)$p.value)
    
    # Skewness 
    normality_test[[paste0(index, "_skewness")]] <- apply(simulated_values[, , index], 1, skewness)
  }
  
  normality_results[[level]] <- normality_test
}

# Stockage final des résultats
normality_results_combined <- bind_rows(normality_results$station, normality_results$depth_layer, .id = "level")

# Combine results
SES_results_depth_layer <- all_results$depth_layer %>%
  tidyr::pivot_longer(cols = -c(name, level), names_to = "index", values_to = "SES_FD") %>%
  mutate(depth_layer=name)

SES_results_station <- all_results$station %>%
  rename(station=name) %>%
  tidyr::pivot_longer(cols = -c(station, level), names_to = "index", values_to = "SES_FD") %>%
  inner_join(metadata %>% select(station, depth), by = "station") %>%
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"
    )
  )

# Combine datasets
SES_results_combined <- bind_rows(SES_results_station, SES_results_depth_layer)

SES_results_combined$depth_layer <- factor(
  SES_results_combined$depth_layer,
  levels = c(
    "Epipelagic",
    "Upper mesopelagic",
    "Lower mesopelagic",
    "Bathypelagic"
  )
)

SES_results_combined$index <- factor(
  SES_results_combined$index,
  levels = c("FRic", "FDis", "FDiv", "FEve"),
  labels = c(
    "Functional richness",
    "Functional dispersion",
    "Functional divergence",
    "Functional evenness"
  )
)

# Plot ----
ggplot(SES_results_combined, aes(x = depth_layer, y = SES_FD)) +
  geom_boxplot(data = SES_results_combined %>% filter(level == "station"),
               aes(color = depth_layer, fill = depth_layer),
               alpha = 0.08, width = 0.5, outlier.shape = NA) +
  geom_point(data = SES_results_combined %>% filter(level == "depth_layer"),
             aes(color = depth_layer),
             size = 5, shape = 18) +
  geom_jitter(data = SES_results_combined %>% filter(level == "station"),
              aes(color = depth_layer),
              size = 1, width = 0.2, alpha = 0.35) +
  geom_hline(yintercept = c(1.96, -1.96), linetype = "dashed", color = "gray40", size = 0.8) +
  scale_color_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
  scale_fill_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
  facet_wrap(~ index, labeller = labeller(index = c(
    FRic = "Functional richness",
    FDis = "Functional dispersion",
    FDiv = "Functional divergence",
    FEve = "Functional evenness"
  ))) +
  labs(x = NULL, y = "Standard Effect Size (SES)") +
  theme_light() +
  theme(
    axis.text.x = element_blank(),
    strip.text.x = element_text(size = 14, color = "black"),
    strip.background = element_rect(fill = "white"),
    axis.title = element_text(size = 13),
    axis.text = element_text(size = 13)
  ) +
  guides(col = "none", fill = "none")

#ggsave("SES_dbFD_two_levels.png", path = "figures", dpi = 700, height = 7, width = 9)

#save(SES_results_combined, normality_results_combined, file = "data.RData")

• Functional dispersion significant in the epipelagic layer

sp_epi <- depth_fish_biomass %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column(var = "depth_layer") %>% 
  filter(depth_layer=="Epipelagic") %>% 
  tidyr::pivot_longer(!depth_layer, values_to = "biomass", names_to = "species") %>% 
  filter(biomass>0) %>% 
  select(species) %>% 
  pull()

alpha_fd_indices_fish <- mFD::alpha.fd.multidim(
  sp_faxes_coord   = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")],
  asb_sp_w         = depth_fish_biomass,
  ind_vect         = c("fdis", "fdiv"),
  scaling          = TRUE,
  check_input      = TRUE,
  details_returned = TRUE)

plots_alpha <- mFD::alpha.multidim.plot(
  output_alpha_fd_multidim = alpha_fd_indices_fish,
  plot_asb_nm              = c("Epipelagic"),
  ind_nm                   = c("fdis", "fdiv"),
  faxes                    = NULL,
  faxes_nm                 = NULL,
  range_faxes              = c(NA, NA),
  plot_sp_nm               = sp_epi,
  save_file                = FALSE,
  check_input              = TRUE) 

plots_alpha$fdis$PC1_PC2

plots_alpha$fdis$PC3_PC4

4. CWM

4.1. CWM booststrap by depth layer

• each trawl is a replica

• non parametric Bootstrap

# Trait 
trait_boot <- morpho_data%>% 
  inner_join(metadata) %>% 
  select(-c(individual_code, years, longitude_start,
            latitude_start, longitude_end, longitude_end,
            volume_filtered, distance_trawled_m)) %>% 
  mutate(
    eye_size = eye_diameter / head_depth,
    orbital_length = eye_diameter / standard_length,
    oral_gape_surface = mouth_width * mouth_depth / body_width * body_depth,
    oral_gape_shape = mouth_depth / mouth_width,
    oral_gape_position = distance_upper_jaws_bottom_head / head_depth,
    lower_jaw_length = lower_jaw_length / standard_length,
    head_length = head_length / standard_length,
    body_depth = body_depth / standard_length,
    pectoral_fin_position = distance_pectoral_bottom_body / body_depth_pectoral_insertion,
    pectoral_fin_insertion = prepectoral_length / standard_length,
    transversal_shape = body_depth / body_width,
    dorsal_fin_insertion = predorsal_length / standard_length,
    eye_position = eye_height / head_depth,
    operculum_volume = operculum_depth / operculum_width,
    gill_outflow = operculum_width,
    caudal_throttle_width = caudal_peduncle_min_depth
  ) %>%
  select(
    depth,
    species,
    eye_size,
    orbital_length,
    gill_outflow,
    oral_gape_surface,
    oral_gape_shape,
    oral_gape_position,
    lower_jaw_length,
    head_length,
    body_depth,
    pectoral_fin_position,
    pectoral_fin_insertion,
    transversal_shape,
    caudal_throttle_width,
    dorsal_fin_insertion,
    eye_position,
    operculum_volume,
    ventral_photophores, 
    gland_head,
    chin_barbel, 
    small_teeth, 
    large_teeth, 
    fang_teeth, 
    retractable_teeth, 
    internal_teeth
  ) %>%
  mutate_at(vars(ventral_photophores, 
                 gland_head,
                 chin_barbel, 
                 small_teeth, 
                 large_teeth, 
                 fang_teeth, 
                 retractable_teeth, 
                 internal_teeth), 
            funs(ifelse(. == "P", 1, ifelse(. == "A", 0, .)))) %>% 
  mutate(across(all_of(c("ventral_photophores", 
                         "gland_head",
                         "chin_barbel", 
                         "small_teeth", 
                         "large_teeth", 
                         "fang_teeth", 
                         "retractable_teeth", 
                         "internal_teeth")), as.numeric)) %>% 
  tidyr::pivot_longer(!c(species,depth), names_to = "trait", values_to = "values")%>%
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))

# Community 
community <-  rbind(data_biomass_2002_2019, data_biomass_2021_2022)%>%
  as.data.frame()%>%
  left_join(metadata) %>%
  select(species, biomass_sp, depth, volume_filtered)%>%
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))%>%
  replace(is.na(.), 0)%>%
  group_by(species, depth)%>%
  mutate(biomass=sum(biomass_sp))%>%
  select(-c(biomass_sp))%>%
  distinct()%>%
  select(-c(volume_filtered)) %>% 
  filter(biomass>0)

trait_filling <- traitstrap::trait_fill(
  # input data (mandatory)
  comm = community,
  traits = trait_boot,
  
  # specifies columns in your data (mandatory)
  abundance_col = "biomass",
  taxon_col = "species",
  trait_col = "trait",
  value_col = "values",
  
  # specifies sampling hierarchy
  scale_hierarchy = c("depth_layer", "depth"),
  
  # min number of samples
  min_n_in_sample = 5
)

# run nonparametric bootstrapping
np_bootstrapped_moments <- traitstrap::trait_np_bootstrap(
  trait_filling, 
  nrep = 100
)

sum_boot_moment <- trait_summarise_boot_moments(
  np_bootstrapped_moments
) %>% 
  mutate(trait= gsub("_"," ", trait)) %>% 
   mutate(trait=stringr::str_to_sentence(trait)) 

sum_boot_moment$depth_layer <- factor(
  sum_boot_moment$depth_layer,
  levels = c(
    "Epipelagic",
    "Upper mesopelagic",
    "Lower mesopelagic",
    "Bathypelagic"
  )
) 

# order traits 
sum_boot_moment$trait <- factor(
  sum_boot_moment$trait,
  levels = c(
    "Caudal throttle width",
    "Oral gape surface",
    "Gill outflow",
    "Large teeth",
    "Eye size",
    "Orbital length",
    "Small teeth",
    "Transversal shape",
    "Body depth",
    "Dorsal fin insertion",
    "Eye position",
    "Oral gape shape",
    "Oral gape position",
    "Internal teeth",
    "Lower jaw length",
    "Pectoral fin position",
    "Ventral photophores",
    "Operculum volume",
    "Pectoral fin insertion",
    "Head length",
    "Chin barbel",
    "Fang teeth",
    "Gland head",
    "Retractable teeth"
  )
) 

ggplot(sum_boot_moment, aes(x = depth_layer, y = mean)) + 
  geom_point(alpha = 0.5, size = 1, position = position_jitter(width = 0.2), aes(col=depth_layer)) +
  geom_boxplot(aes(group=depth_layer, col=depth_layer, fill=depth_layer), alpha=0.1)+
  scale_color_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
  scale_fill_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
  facet_wrap(~trait, scales = "free", ncol = 4) +
  theme_light() +
  labs(y="Community biomass-weighted mean (CWM) ")+
  guides(col="none", fill="none")+
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        panel.grid.minor = element_blank(),
        #panel.grid.major = element_blank(),
        axis.title.y.left = element_text(size = 16),
        axis.title.y = element_text(size = 16),
        axis.text.y = element_text(size = 15),
        strip.background = element_rect(fill = "white"),
        strip.text.x = element_text(size = 16, color = "black"))

ggsave("CWM_boot_sum.png", path = "figures", dpi = 800, height = 13, width = 12)

4.2. PCA CWM

sum_boot_moment_pca <- trait_summarise_boot_moments(
  np_bootstrapped_moments
) %>% 
  ungroup() %>% 
  mutate(trait= gsub("_"," ", trait)) %>% 
  mutate(trait=stringr::str_to_sentence(trait)) %>% 
  select(c(trait, depth_layer, mean)) %>% 
  group_by(
    depth_layer, trait
  ) %>% 
  summarise(median_value= median(mean)) %>% 
  distinct() %>% 
  ungroup() %>% 
  tidyr::pivot_wider(id_cols=depth_layer, values_from = "median_value", names_from = "trait") %>% 
  tibble::column_to_rownames(var = "depth_layer")


res.pca <- FactoMineR::PCA(sum_boot_moment_pca, graph = FALSE)

pca_var <- factoextra::fviz_pca_var(
  res.pca,
  repel = TRUE,
  pointsize = 2,
  arrowsize = 0.1,
  label = "var",
  title = "",
  col.circle = "gray80",
  col.var = "black",
  labelsize=5,
  alpha.arrow=0.5
) +
  labs(x = "", y = "") +
  theme_light() +
  theme(
    aspect.ratio = 1,
    panel.grid.minor =  element_blank(),
    panel.grid.major = element_blank())

pca_ind <- factoextra::fviz_pca_ind(res.pca, label="none", title="") +
  theme_light() +
  xlim(c(-7,7))+
  labs(x="", y="")+
  theme(aspect.ratio = 1)

combined_plot <- gridExtra::grid.arrange(pca_var, pca_ind, ncol = 2)

ggsave("PCA_CWM.png", plot = combined_plot, path = "figures", dpi = 700, height = 6, width = 11)

5. Functional rarity

• Uniqueness - geographical restrectiveness

which traits are mainly driving uniqueness ?

df <- fish_traits %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "species") %>%
  left_join(sp_ui) %>%
  tibble::column_to_rownames(var = "species")

# Identify numerical and categorical traits
numeric_traits <- colnames(df)[sapply(df, is.numeric) &
                                 colnames(df) != "Ui"]

categorical_traits <- colnames(df)[!sapply(df, is.numeric) &
                                     colnames(df) != "Ui"]

# Initialize a data frame to store all results
combined_results_df <- data.frame(
  Trait = character(0),
  Type = character(0),
  Eta_R_squared = numeric(0),
  P_value = numeric(0),
  stringsAsFactors = FALSE
)

# Step 1: Perform Kruskal-Wallis test for categorical traits
for (categorical_trait in categorical_traits) {
  kruskal_result <- kruskal.test(df$Ui ~ df[[categorical_trait]])
  
  # Calculate eta-squared for Kruskal-Wallis
  n_groups <- length(unique(df[[categorical_trait]]))
  n_total <- length(df$Ui)
  h_value <- kruskal_result$statistic
  eta_squared <- (h_value - (n_groups - 1)) / (n_total - n_groups)
  
  combined_results_df <- rbind(
    combined_results_df,
    data.frame(
      Trait = categorical_trait,
      Type = "Categorical",
      Eta_R_squared = as.numeric(format(eta_squared, scientific = FALSE)),
      P_value = as.numeric(format(kruskal_result$p.value, scientific = FALSE)),
      stringsAsFactors = FALSE
    )
  )
}

# Step 2: Fit linear models for numerical traits
for (numeric_trait in numeric_traits) {
  lm_result <- lm(Ui ~ df[[numeric_trait]], data = df)
  summary_stats <- summary(lm_result)
  
  combined_results_df <- rbind(
    combined_results_df,
    data.frame(
      Trait = numeric_trait,
      Type = "Numeric",
      Eta_R_squared = as.numeric(format(summary_stats$r.squared, scientific = FALSE)),
      P_value = as.numeric(format(
        summary_stats$coefficients[2, 4], scientific = FALSE
      )),
      stringsAsFactors = FALSE
    )
  )
}

# Round the numeric columns to three decimal places
combined_results_df[, c("Eta_R_squared", "P_value")] <-
  round(combined_results_df[, c("Eta_R_squared", "P_value")], 3)

htmltools::tagList(DT::datatable(combined_results_df))

6. Appendices

# Convert sp_faxes_coord_fish to a data frame
sp_coord_community_2 <- as.data.frame(sp_faxes_coord_fish[, c("PC3", "PC4")]) %>%
  tibble::rownames_to_column(var = "species")

# Create a list to store data frames for each depth layer
all_data_layers_2 <- lapply(unique(sp_all_layers$depth_layer), function(layer) {
  sp_layer <- sp_all_layers %>%
    filter(depth_layer == layer)
  
  #presence or absence
  sp_coord_layer <- sp_coord_community_2 %>%
    mutate(layer_presence = ifelse(species %in% sp_layer$species, "yes", "no")) %>%
    mutate(depth_layer = layer)
  
  return(sp_coord_layer)
})

# Combine all layers into one data frame
sp_coord_community_all_layers <- bind_rows(all_data_layers_2)

# Function to calculate the convex hull for a given data frame
calculate_hull <- function(data) {
  if (nrow(data) < 3) return(data)  
  data %>%
    slice(chull(PC3, PC4))
}

# Initialize lists to store hulls
hull_all_combined <- NULL
hull_layer_combined <- NULL

for(layer in unique(sp_coord_community_all_layers$depth_layer)) {
  sp_layer <- sp_coord_community_all_layers %>%
    filter(depth_layer == layer)
  
  hull_all <- calculate_hull(sp_layer) %>%
    mutate(depth_layer = layer, type = "all")
  
  hull_layer <- calculate_hull(sp_layer %>% filter(layer_presence == "yes")) %>%
    mutate(depth_layer = layer, type = "present")
  
  hull_all_combined <- bind_rows(hull_all_combined, hull_all)
  hull_layer_combined <- bind_rows(hull_layer_combined, hull_layer)
}

# Define depth layer levels
depth_levels <- c("Epipelagic", "Upper mesopelagic", "Lower mesopelagic", "Bathypelagic")

sp_coord_community_all_layers$depth_layer <- factor(sp_coord_community_all_layers$depth_layer, levels = depth_levels)
hull_all_combined$depth_layer <- factor(hull_all_combined$depth_layer, levels = depth_levels)
hull_layer_combined$depth_layer <- factor(hull_layer_combined$depth_layer, levels = depth_levels)

# Define colors for each depth layer
depth_colors <- c("Epipelagic" = "#FEA520", 
                  "Upper mesopelagic" = "#D62246", 
                  "Lower mesopelagic" = "#6255B4", 
                  "Bathypelagic" = "#3C685A")

# Plot
plot_depth <- ggplot() +
  # Plot present species
  geom_point(data = sp_coord_community_all_layers %>% filter(layer_presence == "yes"), 
             aes(PC3, PC4, col = depth_layer), size = 2, shape=19) +
  # Plot absent species
  geom_point(data = sp_coord_community_all_layers %>% filter(layer_presence == "no"), 
             aes(PC3, PC4), col = "grey", shape = 8, size = 2) +
  # Plot hulls
  geom_polygon(data = hull_all_combined, aes(x = PC3, y = PC4), 
               fill = "gray70", alpha = 0.2) +
  geom_polygon(data = hull_layer_combined, aes(x = PC3, y = PC4, group = depth_layer, fill = depth_layer), 
               alpha = 0.4) +
  theme_light() +
  scale_color_manual(values = depth_colors) +
  scale_fill_manual(values = depth_colors) +
  guides(col = "none", shape = "none", fill = "none") +
  facet_wrap(~ depth_layer, nrow=1 ) +  
  theme(
    strip.text.x = element_text(size = 14, color = "gray20"),
    strip.background = element_rect(fill = "white"),
    aspect.ratio = 1,
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank()
  )


beta_fd_indices_fish <- mFD::beta.fd.multidim(
  sp_faxes_coord   = sp_faxes_coord_fish[ , c("PC3", "PC4", "PC3", "PC4")],
  asb_sp_occ       = asb_sp_fish_occ,
  check_input      = TRUE,
  beta_family      = c("Jaccard"),
  details_returned = T)

Serial computing of convex hulls shaping assemblages with conv1

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
Serial computing of convex hulls shaping assemblages with conv2

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
Serial computing of intersections between pairs of assemblages with inter_geom_coord

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |======================================================================| 100%
Serial computing of intersections between pairs of assemblages with inter_rcdd_coord & qhull.opt1

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |======================================================================| 100%
Serial computing of intersections between pairs of assemblages with inter_rcdd_coord & qhull.opt2

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |============                                                          |  17%

vertices_community <-beta_fd_indices_fish$"details"$"asb_vertices"$Bathypelagic 

sp_coord_community_2 <- as.data.frame(sp_faxes_coord_fish[, 3:4]) %>% 
  tibble::rownames_to_column(var = "species") %>% 
  mutate(vertices=case_when(species%in%vertices_community~"vertice",
                            !species%in%vertices_community~"not_vertice")) 

hull <- sp_coord_community_2%>% 
  slice(chull(PC3, PC4))

plot_com <- ggplot(data = sp_coord_community_2,aes(PC3, PC4)) +
  scale_color_manual(values = c("grey", "black"))+
  geom_point(size = 2.5,aes( shape=vertices, col= vertices)) +
  geom_polygon(data=hull, aes(x = PC3, y = PC4), 
               fill = "gray70", alpha = 0.2) +
  theme_light() +
  geom_hline(yintercept = 0, linetype = 2, col = "gray45") +
  geom_vline(xintercept = 0, linetype = 2, col = "gray45") +
  guides(col = "none", shape = "none", fill = "none") +
  theme(
    strip.text.x = element_text(size = 14, color = "gray20"),
    strip.background = element_rect(fill = "white"),
    aspect.ratio = 1,
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank()
  )

ggpubr::ggarrange(
  labels = c("A","B"),
  plot_com,   
  plot_depth, 
  nrow = 2
) 

ggsave("functional_space_PC3_4.png", path = "figures", dpi = 800, height = 8, width = 10)